Theory of spin-phonon coupling in multiferroic Mn perovskites 
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Magnetoelectric phase diagrams of the rare-earth (R) Mn perovskites J?Mn03 are theoretically 
studied by focusing on crucial roles of the symmetric magnetostriction or the Peierls-type spin- 
phonon coupling through extending our previous work [M. Mochizuki et al., Phys. Rev. Lett. 
105, 037205 (2010)]. We first construct a microscopic classical Heisenberg model for i?MnC>3 in- 
cluding the frustrated spin exchanges, single-ion anisotropy, and Dzyaloshinskii-Moriya interaction. 
We also incorporate the lattice degree of freedom coupled to the Mn spins via the Peierls-type 
magnetostriction. By analyzing this model using the replica-exchange Monte-Carlo technique, we 
reproduce the entire phase diagram of i?Mn03 in the plane of temperature and magnitude of the 
orthorhombic lattice distortion. Surprisingly it is found that in the a&-plane spiral spin phase, the 
(S ■ S)-type magnetostriction plays an important role for the ferroelectric order with polarization 
P\\a whose contribution is comparable to or larger than the contribution from the (S x S')-type 
magnetostriction, whereas in the bc-plane spiral phase, the ferroelectric order with P\\c is purely of 
(S x S) origin. This explains much larger P in the afc-plane spiral phase than the fcc-plane spiral 
phase as observed experimentally, and gives a clue how to enhance the magnetoelectric coupling in 
the spin-spiral-based multiferroics. We also predict a noncollinear deformation of the _B-type spin 
structure resulting in the finite (S x S) contribution to the ferroelectric order with P\\a, and a wide 
coexisting regime of the commensurate E and incommensurate spiral states, which resolve several 
experimental puzzles. 

PACS numbers: 75.80.+q, 75.85.+t, 75.47.Lx, 75.10.Hk 



I. INTRODUCTION 

Since the magnetic (electric) induction of electric po- 
larization (magnetization) was proposed theoretically by 
Dzyaloshinskii in 1959 1 -, the magnetoelectric coupling in 
solids has attracted a great deal of interest. Many mag- 
netic materials have been demonstrated to exhibit the 
magnetoelectric effect^, but the observed effect was very 
weak. Recently the interest has been revived by discov- 
ery of the magnetically induced ferroelectric order, i.e., 
multiferroic order in perovskite TbMnOa 4 . 

An innovative aspect of this discovery is that in 
TbMnC>3, although the lattice structure retains the in- 
version symmetry, a nontrivial magnetic order breaks 
the inversion symmetry and induces the ferroelectric po- 
larization^—. This is in striking contrast to the usual 
ferroelectrics whose ferroelectricity originates from the 
crystal structure with inherent broken inversion sym- 
metry. Therefore the magnetoelectric coupling is very 
strong in TbMnC>3, which leads to a lot of intriguing 
cross-correlation phenomena such as electromagnon ex- 
citations^—, magnetic-field control of the ferroelectric- 
it\— — — , colossal magnetocapacitanc o 18 ' 23 " — and so on. 

TbMnC>3 shows successive two magnetic phase transi- 
tions with lowering temperature, and below the second 
transition, the ferroelectric polarization P appears and 
grows as temperature decreases. The emergence of P in 
this compound is explained by the antisymmetric mag- 
netostriction associated with the cross-product of spins 



(Si x Sj) as described by2£2£ 
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A <■>.> x {Si xSj). 



(1) 
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Here e^j is the unit vector connecting two spin sites i and 
j, and A is a coupling constant determined by the spin- 
exchange and spin-orbit interactions. This formula im- 
plies that two canted spins Si and Sj can induce an elec- 
tric polarization pij via the spin-orbit coupling. Conse- 
quently a transverse spiral spin order as a sequence of the 
canted spins can generate ferroelectric Pas- A neutron- 
scattering experiment for TbMnC>3 confirmed that the 
Mn spins in its multiferroic phase with P\\c rotate within 
the be plane to form a transverse spin spiral propagating 
along the b axis 2 ^. For the a, b and c axes, we adopt the 
Pbnm setting [see Fig. |3J^a)]. 

This equation also implies that the direction of Pas 
depends on orientation of the spin spiral plane. In the 
frc-plane (a6-plane) spiral spin order, the Pas directs 
in the c (a) direction as shown in Fig. [TJ This rela- 
tionship has been confirmed by neutron-scattering ex- 
periments 2 ^—. In i?Mn03 with R being a rare-earth 
ion, the spiral-plane orientation is determined by a sub- 
tle competition between the magnetic anisotropy and 
the Dzyaloshinskii-Moriya interactio n 33 ' 34 . The ground 
states of TbMnC>3 and DyMnC>3 exhibit the 6c-planc 
spiral order with Pas||c, while the a6-plane spiral or- 
der with Pas || a arc observed in some solid solutions 
Eui^Y^MnOs and Gdi-^Tb^MnOs. 
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FIG. 1: (Color online) (a) Relationship between the spiral- 
plane orientation and the spontaneous ferroelectric polariza- 
tion Pas in the foe-plane spiral spin structure predicted by the 
spin-current model-^"— ^. (b) That in the afo-plane spiral spin 
structure. 



Another multifcrroic phase was theoretically pre- 
dicted^— and experimentally discovered 3 ^— in 
.RM11O3 with much smaller R ions such as Y, Ho, Tm, 
Lu as well. These compounds exhibit the E-type 
antifcrromagnctic ground state where the Mn spins 
form an up-up-down-down structure. In this state, 
the symmetric magnetostriction associated with the 
inner-product of the spins (Si ■ Sj) induces a ferroelectric 
polarization Pg parallel to the a axi o 31 ' 42 , which is given 
by 



<ij> 



& -St). 



(2) 



Here TZij is a form factor which reflects the zigzag MnO 
chains, and is nonzero because of the absence of inver- 
sion symmetry at the center of Mn-O-Mn bond. It is 
worth mentioning that Pg in the E-typc phase is much 
larger in magnitude than PAg in the spiral spin phases be- 
cause the symmetric-magnetostriction mechanism is as- 
sociated only with the spin-exchange interaction J but 
not with the spin-orbit interaction or the Dzyaloshinskii- 
Moriya interaction in contrast to the antisymmetric- 
magnetostriction mechanism. Typically the magnitude 
of Pg in the E-type phase is ~4600 /iC/m 2 , whereas that 
of PAg in the spiral spin phases takes ^500 /iC/m 2 at 
most in i?MnC>3. 

The perovskite structure of .RM11O3 is orthorhombi- 
cally distorted with alternately tilted MnOg octahedra. 
Magnitude of this GdFe03-type distortion varies depend- 
ing on the size of the R ion. With a smaller R ion, the 
lattice is more significantly distorted, and the Mn-O-Mn 
bond angle is reduced more from 180°. The magnetoelec- 
tric phase diagram of .RM11O3 as a function of the magni- 
tude of the GdFc03-type distortion or the ionic i?-site ra- 
dius (vr) has been studied experimental ^ 38 ' 43 " —. It has 
been revealed that following four magnetoelectric phases 
successively emerge at low temperatures with decreasing 
rpr^ [see Figs.^a) and (b)]; (i) non-ferroelectric ^4-type 
phase where the Mn spins align fcrromagnctically in the 
ab plane, (ii) ferroelectric a6-plane spiral spin phase with 
P||a, (iii) ferroelectric 6c-plane spiral spin phase with 
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FIG. 2: (Color online) (a) Experimentally obtained mag- 
netoelectric phase diagram of _RMn03 and solid solutions 
Gdi-^Tb^MnOs and (b) that of solid-solution systems 
Eui-^YzMnOs and Yi-^LuyMnOs in the plane of tempera- 
ture and (effective) ionic radius of the R ion^. 



P||c, and (iv) ferroelectric -E-type phase with very large 
P\\a. In all these phases, the Mn spins are stacked in 
a (nearly) staggered manner along the c axis because of 
the strong interplane antiferromagnetic coupling. 

On top of these four phases, there exists a paraelec- 
tric sinusoidal collincar spin phase in the intermediate 
temperature regime where the collinear Mn spins paral- 
lel to the b axis are sinusoidally modulated in amplitude. 
A theoretical model for i?Mn03 which thoroughly de- 
scribes severe competitions among these magnetoelectric 
phases has long been desired. 

Here we briefly introduce previous theoretical stud- 
ies on this issue. In the canonical two-dimensional 
•h-J-2 classical Heisenberg model with ferromagnetic 
nearest-neighbor interaction Ji(<0) and antiferromag- 
netic second-neighbor interaction </2(>0) (see Fig. 18 
of Reff 3 ^), the ferromagnetic order is realized for 
|</2/</i|<0.5, while the spiral order is stabilized for 
l^/ Ji I >0.5. In i?Mn03, the second-neighbor exchange 
J2 originates from an indirect overlap of the Mn 3c? or- 
bitals via the inbetween two O 2p orbitals enhanced by 
the orthorhombic lattice distortion, and thus is de- 
pendent^. This can explain the observed phase evolu- 
tion from the ^4-type phase (i.e., staggered stacking of the 
ferromagnetic planes) to the spiral phase with decreasing 
t\r. Within this model, however, following experimental 
observations in i?Mn03 cannot be reproduced: (i) Stabil- 
ities of the specific spin-spiral planes (ab- or frc-cycloidal 
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spin structures), which depend on r^, temperature, and 
magnetic field, (ii) Emergence of the E-type phase ac- 
companied by ferroelectricity. 

In order to reproduce the parameter-dependent spiral- 
plane directions, magnetic anisotropics are essentially im- 
portant. Taking into account the orthorhombic lattice 
structure of i?Mn03, it is quite natural to examine the 
single-ion anisotropy and the Dzyaloshinskii-Moriya in- 
teraction reflecting the local lattice distortion. By in- 
troducing such interactions to the frustrated classical 
Heisenbcrg mode l 22 ' 33 ' 34 , the authors have successfully 
reproduced the tr-T and T-H phase diagrams of i?Mn03 
with respect to the spin structures. It is noteworthy that 
the observed 90° flop of the spin-spiral plane between 
ab and be has been ascribed to the competition between 
the single-ion anisotropy and the Dzyaloshinskii-Moriya 
interaction, and thus both ingredients are indispensable. 

This model, however, fails to reproduce the E-type 
phase. An attempt to understand the E-type phase was 
performed using an Ising model 4 ^, which however failed 
to explain the spiral phases by definition. A two-orbital 
double-exchange model was examined, and the E-type 
spin order as well as the phase evolution from spiral to E- 
type phases were reproduced 4 ^. In the double-exchange 
model, however, the electron correlation in the i2g-orbital 
sector is neglected although it is usually strong enough to 
make the system Mott insulating. Thus the mechanism of 
the E-type order in the double-exchange model may not 
be straightforwardly applicable to the present undoped 
i?Mn03 system. Kaplan et al proposed a biquadratic 
spin interaction originating from the spin-phonon cou- 
pling as an origin of the E-type orde r 50 ' 51 . Bond al- 
ternation or staggered modulation of the ferromagnetic 
exchanges was also proposed for its origin 5 ^. The latter 
two works suggest the importance of the spin-lattice cou- 
pling to understand the emergence of E-typc order and 
its ferroelectricity. After these attempts, the entire phase 
diagram of i?MnC>3 including all the competing magneto- 
electric phases was reproduced by a classical Heisenberg 
model including the Peierls-type spin-lattice coupling 5 ^. 
This spin-lattice coupling is a source of the ferroelectric 
polarization associated with the symmetric magnetostric- 
tion given by Eq. $Z§. 

Here we note that the biquadratic interaction, 
— Ebiq Yl<i j>(&i ■ Sj) 2 , is an effective interaction among 
spins via the spin-phonon coupling, which is derived by 
integrating out the phonon degrees of freedom. This in- 
teraction favors the collinear spin alignment, and thus 
stabilizes the E-type spin order as compared to the cy- 
cloidal orders. However, since this interaction does not 
contain the phonons or the lattice degrees of freedom ex- 
plicitly, it is not appropriate to study behaviors of the 
lattice displacements or the electric polarizations as well 
as the multiferroic properties. In contrast, the Peierls- 
type spin-phonon coupling is adequate for such studies. 

In this paper, we study theoretically origins and prop- 
erties of the magnctoclcctric phases in i?MnC>3 by focus- 
ing on roles of the Peierls-type spin-phonon coupling on 



the basis of the Monte-Carlo analysis of a spin model. 
We reveal a large contribution of the (S ■ S')-type magne- 
tostriction to P\\a in the afo-plane spiral phase in addi- 
tion to the (S x 5')-type magnetostriction. This finding is 
quite surprising because the (S x S) mechanism has been 
considered as a unique origin of the ferroelectric polariza- 
tion in the spin spiral phase thus far. On the other hand, 
the P\\c in the &c-plane spiral phase is purely of (S x S) 
origin. This solves a long standing puzzle of much larger 
P observed in the a6-plane spiral phase. This (S ■ S) 
mechanism can be generally expected in many other spin- 
spiral-based multiferroics, and gives a clue how to design 
the enhanced magnetoelectric coupling in materials. We 
also predict a cycloidal deformation of the E-type spin 
structure, which causes an additional (S x S) contribu- 
tion to the ferroelectric order with -P||a, in addition to 
the dominant (S ■ S) contribution. In addition, we find a 
wide regime where the E-type and spiral states coexist. 
On the basis of these findings, we resolve a puzzle in the 
neutron-scattering experiments for i?Mn03 with i?=Y, 
Ho, and Er. 

The rest of this paper is organized as follows. In Sec. II, 
we introduce the classical spin model for /2M11O3 includ- 
ing the Peierls-type spin-phonon coupling. Then we ex- 
plain the methods for numerical simulations, and phys- 
ical quantities which we calculate in the simulations in 
Sec. III. In Sec. IV, we discuss the results for the whole 
phase diagram, the spiral spin states, the E-type state, 
and the coexistence of the E-type and incommensurate 
states in each subsection. Section V is devoted to the 
summary. A short report of the present work has been 
published 5 ^. In addition to the detailed explanation, 
some further results are presented in this paper. 



II. MODEL 

To describe the Mn 3d-spin system in i?Mn03, we em- 
ploy a classical Heisenberg model on a cubic lattice 5 ^, in 
which the Mn S —2 spins arc treated as classical vectors, 
S^i^/S 2 -S'i cos 9 t , ^S 2 -S 2 t sm6 u S cl ) with respect 
to the a, b, and c axes. The Hamiltonian is given by 

T~L = Hex + 7~t-Ea + ^ifa + ^dm + Hk, (3) 

with 
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— JijSi • Sj, 

<i,j> 


(4) 
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FIG. 3: (Color online) (a) Spin-exchange interactions in 7?MnC>3 and tilted local coordinate axes £j, r/j, and C,i attached to the 
ith MnOe octahedron. Here FM (AFM) denotes (anti)ferromagnetic exchange interaction. For the spin-exchange interactions, 
we consider ferromagnetic exchange J a t on the Mn-Mn bonds along the pseudocubic x and y axes, (anti)ferromagnetic exchange 
Ja {Jb) on the in-plane diagonal Mn-Mn bonds along the a (b) axis, and antiferromagnetic exchange J c along the c axis, (b) 
Mn(i)-0-Mn(j) bond in the orthorhombic lattice and local vector riij. The O ion is displaced from its cubic position (0) to the 
orthorhombic position (A ) at higher temperatures. At low temperatures, a further shift 5ij along riij can be induced by the 
spin-lattice coupling in the presence of magnetic order, (c) A vs J a t for several i?Mn03 compounds calculated in Ref.— , which 
gives J' a i=d Jab /9A =2.5. Here A is normalized by the MnO bond length, (d) Main exchange path for the next-neighbor 
ferromagnetic exchange J a along the a axis [(blue) dotted line]. 



TABLE I: Structural parameters of DyMnOs from Refi* 



(A) b (A) c {A) XOl y Ql 



XQ 2 



yo 2 



zo 2 



5.2785 5.8337 7.3778 0.1092 0.4642 0.7028 0.3276 0.0521 



where i Xl i y and i z represent the integer coordinates of 
the ith Mn ion with respect to the pseudocubic x, y and 
z axes [see Fig.[3ja)]. 

The first term H ex describes the spin-exchange inter- 
actions as shown in Fig. EJa). Since the strength of the 
nearest-neighbor ferromagnetic coupling in i?MnC>3 sen- 
sitively depends on the Mn-O-Mn bond angle, we con- 
sider the Peierls-type spin-phonon coupling 



Jij — Jab 



Jab&ijj 



(9) 



for the in-plane Mn-O-Mn bonds where J' ab =dJ a b/ 98. 
Here Sij and 6 denote a shift of the O ion between 
ith and jth Mn ions normalized by the averaged MnO 
bond length. Note that the O ion in the orthorhom- 
bic lattice is already displaced from its cubic position. 
We consider 5ij as a further shift of the position in the 
presence of magnetic order at low temperatures with re- 
spect to the orthorhombic position at higher tempera- 
tures. We assume that the shift of the O ion Sij oc- 
curs along the local axis n.ij directing from its cubic po- 
sition (0) to the orthorhombic position (A Q ) at higher 
temperature as shown in Fig. EJb). Then the positive 
(negative) shift decreases (increases) the Mn-O-Mn bond 
angle. The nearest neighbor ferromagnetic exchange be- 
comes stronger (weaker) as the Mn-O-Mn bond angle in- 
creases (decreases), which implies positive J' ab . 



The second and the third terms, H5„ and TLS a , stand 

' old aid ' 

for the single- ion anisotropics. Here 77.; and Q are 
tilted local axes attached to the ith MnOe octahedron 
as shown Fig. [3ja). The former term makes the mag- 
netization along the c axis hard, while the latter term 
causes alternation of the local easy and hard magnetiza- 
tion axes along the £j and rji axes in the ab plane owing 
to the staggered d^ x 2_ r i. / d^ y 2_ r 2 type orbital ordering. 
The directional vectors r/j and with respect to the 
a, b and c axes are given by 



& = 



a[0.25 + (-iy*+ l « (0.75 - xo 2 )] 
&[0.25-(-l) < -+**(y Oa -0.25)] 

a[-0.25+(-l)^ +i "(0.75-a;o 2 )] 
6[0.25+(-l)^+Myo 2 -0.25)] 



-c(-l) 



ZQ 2 



6(-l)^(0.5-2/ Ol ) 
0.25c 



(10) 



(11) 



(12) 



Here jo 2 , yo 2 an d z o 2 ( x Oi an d UOi) are the coordina- 
tion parameters of the in-plane (out-of-plane) oxygens, 
and a, b and c are the lattice parameters. For values 
of these parameters, we use the experimental data of 
DyMnO^ throughout the calculations (see Table [T|. 

The fourth term, Hdm, denotes the Dzyaloshinskii- 
Moriya interaction^—. The vectors are defined on 
the Mn(i)-0-Mii(j) bonds. Because of the crystal sym- 
metry, they are expressed using five parameters, a a b, /3 a b, 
lab, etc-, and /3 C , as given in Refi^. Their expressions are 
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drriT = (aab, Rab, yab) 




.dm5 = (-aab, Rab, yab) 



I : °5 



dmS = (aab, Raft, yab) V 
o 

dl2 = (aab, -Rab, yab) 



or" 

2 



dim = (-ac, -Rc, 0) 
Q^- d/3 = (-aab, -Rab, yab) 



r: d/1 = (aab, -Rab, yab) 

dl4 = {-aab, -Rab, yab) ^ Q 



FIG. 4: (Color online) Dzyaloshinskii-Moriya vectors dij as- 
sociated with the Mn(i)-0-Mn(j) bonds. 



TABLE II: Model parameters used in the calculations. The 
energy unit is meV. 



rlex 


Jab = -0.8, Ja = -0.1, 




1.25 J^ = 2.5 


niD njE 
''•siai 'T-sia 


D=0.2, 


£■=0.25 






?<DM 


a ab =0.1, 


/3ab=0.1, 


7ab = 


=0.14 




a c =0.42, 


/3 C =0.1 






Hk 


K=500 









[see Fig. [3jc)] , which gives J \ b =dJ a b/dA =2.5 meV. 

The main exchange path for the ferromagnetic J a is 
shown in Fig. EJd), which contains two O 2p orbitals. 



The electron starting from the Mnl d 



% 2 



orbital dom- 



inantly reaches the unoccupied Mn2 d z 2_ x 2 orbital via 
the 2p x and 2p y orbitals. This results in the doubly occu- 
pied Mn2 ion with d^ y2 _ r 2d 1 z2 _ x 2 electron configuration 
as an intermediate state of the perturbation process with 
respect to the d-p and p-p transfer integrals. This pro- 
cess favors the parallel spin configuration relative to the 
antiparallel one since the Hund's-rule coupling reduces 
energy of the intermediate state, which leads to the fer- 
romagnetic exchange J a . 

We treat the antiferromagnetic exchange J b as a vari- 
able which increases (decreases) as tr decreases (in- 
creases). This is because the exchange path for Jb con- 
tains two O 2p orbitals between two Mn e g orbitals neigh- 
boring along the b direction, and the orthorhombic dis- 
tortion, whose magnitude is controlled by r^, enhances 
their p-p hybridizatio n 48 ' 59 . On the other hand, previous 
microscopic evaluations of the model parameters showed 
that the parameters except for Jb are almost insensitive 
to the i?-site species in/near the multifcrroic phases^. 
We find that overall features of the phase evolution upon 
the .R-site variation can be reproduced as a function of 
Jb even without considering the slight independence of 
other parameters. 



pven by (see also Fig. 0]), 



III. METHOD 
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*Pab 
lab 

OLab 
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(_2_)'ix+i H 
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i,i-\-z 



(-1) 

(-1)*-+**+*- /?„(, 





(13) 



(14) 



(15) 



The last term represents the lattice elastic term with K 
being the elastic constant. 

The values of J a b, J c , Jb, D, E, and five Dzyaloshinskii- 
Moriya parameters have been microscopically determined 
or have been estimated in Ref^ for several _RMn03 com- 
pounds. Except for J b , they are nearly invariant upon the 
i?-site variation in the vicinity of the multifcrroic phases. 
We also find that very weak ferromagnetic exchange J a 
along the a axis is necessary to produce the -E-type phase, 
and adopt J a =— 0.1 meV. The value of K is chosen so as 
to reproduce the experimental value of P (~4600 pC/m 2 ) 
in the -E-type phase, which mostly comes from the (S ■ S) 
contribution. The model parameters used in the calcula- 
tions are summarized in Table [TT1 We obtain the value of 
J' ab from the A D dependence of J ab for several R species 



We investigate finite-temperature properties of the 
Hamiltonian by using the Monte-Carlo technique. 
The spin system in i?Mn03 has frustrating interactions 
and exhibits first-order phase transitions, so that it is 
difficult to treat this system using conventional serial- 
temperature Monte-Carlo methods. Thus we adopt the 
replica-exchange Monte-Carlo method^, in which one 
simulates N-r replicas at different temperatures cover- 
ing not only a low-temperature regime of interest but 
also a higher-temperature regime above it, and allows 
configurational exchange between the replicas. The in- 
clusion of high-temperature configurations enables the 
lower-temperature systems to access a broad phase space 
and to avoid being trapped in local energy minima. We 
perform the simulations for temperature range 0.5<fcBT 
(meV)<7 with 200 temperature meshes. Both spins and 
oxygen positions (Sij) are updated and relaxed in the 
simulation. To achieve efficient updates of the oxygen 
displacements, we adopt a window for Monte-Carlo sam- 
pling of Si j, —W < St j < W with W=0.5, within which 
we generate random numbers for the Sij sampling. We 
carry out each configurational exchange after every 400 
standard Monte-Carlo steps. Typically, we perform 1000 
exchanges after the sufficient thcrmalization steps for sys- 
tems with 7V=48x48x6 sites along the x, y and z axes 
under the periodic boundary condition. 

We identify transition points and magnetic structures 
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by calculating temperature profiles of the specific heat 
C S (T) for the spin system and the 7-axis component of 
the total spin-helicity vector h 7 (T) with 7=a, b and c. 
They are respectively calculated by 

C S (T) = ±d(H-H K )/d(k B T), (16) 
K(T) = ^(| ^(S, x S t+£ + S t x S l+ y) 7 \)/S 2 . 

i 

(17) 

Here the brackets denote the thermal average. 

Using the point-charge model, we calculate the electric 
polarization Ps = (P a , Pf,, P c ) due to the oxygen displace- 
ments 5ij induced by the (S ■ S^-type magnetostriction. 
Considering the staggered arrangement of the local axes 
riij on the the zigzag Mn-0 chain, we calculate E 7 ("f=a, 
b, c) by 

p i = -^(iE[(- 1 ) Jl+i!,+Tn ^^+(- 1 ) la!+l! ' +,l(5 M +2 ,]i} ; 

(18) 

where (to, n)=(0, 0) for 7=a, (to, n)=(l, 0) for 7=6, 
and (to, n)=(i 2 +l, i z +l) for 7=c. Here the constant n 7 
is calculated to be 4.6 xlO 5 /zC/to 2 for j—a and b, and 
4.7xl0 5 /j,C/m 2 for 7=0 from the lattice parameters. 

In order to confirm the magnetic structures, we also 
calculate the spin and spin-helicity correlation functions 
in the momentum space, S* 7 (fc, T) and H 7 (k, T) for 7=0, 
&, and c. They are calculated by 

S 7 (fc,T) = -l^Y.iSr^e^^K (19) 

id 
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FIG. 5: (Color online) Theoretical phase diagram of i?Mn03 
in the plane of temperature and Jb- Here ICS denotes the 
incommensurate spiral phase. In the shaded area, incommen- 
surate spin states can coexist with the .E-type state (see text). 



Here ac is the dimcnsionlcss Gilbert-damping coefficient 
introduced phcnomenologically. For the value of ac, we 
take a rather small value of «g = 0.01 to achieve a 
slow relaxation towards a real ground-state spin structure 
with a minimum energy. We solve this equation using 
the fourth-order Runge-Kutta method after the linear- 
ization. For the convergence, we use thermally relaxed 
spin configurations obtained in the Monte-Carlo simula- 
tions at low temperatures as initial states. 



IV. RESULTS 
A. Phase Diagram 



1 -j 



Here h 7 i is the 7 component of the local spin-helicity 
vector hi — (h a i, h b i, h C i), which is defined as 

hi = ^(Si x S l+£ + S t x S l+ y)/S 2 . (21) 

In the following, we write these correlation functions sim- 
ply as Sj(k) and H 1 (k) by omitting T. 

On the other hand, we study ground-state properties of 
the Hamiltonian ([3]) by numerically solving the Landau- 
Lifshitz-Gilbert equation; 



dt ~ blX + s SlX dt 



The effective local magnetic fields H° s acting on the ith 
Mn spin Si are derived from the spin-derivative of the 
Hamiltonian % as 



Hf = -dU/dSi. 



In Fig. [5j we display theoretically obtained phase di- 
(r-i-r-j) ^0) agram in the plane of temperature k^T and antiferro- 
magnetic exchange along the b axis Jy^, which suc- 
cessfully reproduces the experimental phase diagram in 
Figs. [2ja) and (b). At low temperatures, the A-type, 
a&-plane spiral, fee-plane spiral, and -E-type phases suc- 
cessively emerge as J b increases. Here the magnetic 
structure is commensurate with g;,=0.57r in the E-type 
phase, whereas it is incommensurate in the ab and bc- 
plane spiral phases. The sinusoidal collinear state is also 
incommensurate even above the E-type phase, and the 
spin-phonon coupling is a source of the incommensurate- 
commensurate transition with lowering temperature. In 
the shaded area, although the E-type state has the low- 
est energy, incommensurate spin states have deep energy 
minimum, and can coexist with the E-type state as will 

(22) be discussed in Sec.IV-D. 
Note that the experimental phase diagram of 

the solid-solution systems, i.e., Eui-^Y^MnOs and 
Yi_j,Lu. y Mn03 in Fig. [2jb) shows that on the verge of 
the phase boundary between the 6c-plane spiral and the 
E-typc phases, the transition temperature is strongly 

(23) suppressed, and a V-shaped bicritical point appears. 
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FIG. 6: (Color online) Temperature profiles of specific heat 
C S (T), spin helicity /i 7 (T), and polarizations due to the (S ■ 
S)-type magnetostriction P-,(T) for J{,=2.4 meV. 



This can be attributed to the randomness effect inher- 
ent to the solid-solution systems, which suppresses the 
transition temperature of the first-order phase transi- 
tion. Contrastingly, in the experimental phase diagram 
of the compounds RMnOs with almost no randomness 
effect [Fig. HJa)], only a very small dip appears. The 
randomness-induced suppression of the first-order tran- 
sition temperature at the bicritical point has been es- 
tablished well by the early experimental and theoretical 
studies on the solid-solution Mn pcrovskite o 61 ' 62 . On the 
other hand, in the theoretical phase diagram in Fig. [SJ 
which is obtained without considering the randomness 
effect, shows a nearly straight phase boundary with no 
anomaly. 



B. i5-type Spin Phase 

First we discuss the E-type spin phase. In Fig. |B1 
we display calculated temperature profiles of the specific 
heat C s (T) , the spin helicity hy (T) , and the polarization 
due to the (S ■ S)-type magnetostriction P 7 (T) for Jf,=2.4 
meV. The system exhibits two phase transitions and 
three magnetic phases emerging successively with low- 
ering temperature, i.e., the paramagnetic, the sinusoidal 
collinear, and the B-type phases^. The magnetic struc- 
ture is incommensurate (qb=0.4587r) in the sinusoidal 
collinear phase, while it is commensurate (qb—0.5ir) in 
the B-type phase. 

Interestingly we find a finite c-axis component h c (T) 
of the spin helicity in the B-type phase, indicating that 
its spin structure is not collinear in reality, but its up- 
up-down-down structure is subject to a cycloidal defor- 
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FIG. 7: (Color online) (a) Real-space spin configuration of the 
-B-type phase for two kinds of ab planes, Zi~0 and 0.5. (b) 
Spin alignment on the zigzag Mn-0 chain along the x axis. 
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FIG. 8: (Color online) (a) Pure collinear up-up-down-down 
spin structure with Mn spins parallel to the b axis, (b) 
Because of the staggered d- i( 2 _ r n / 'd 3r .2 _ r 2 orbitals, the lo- 
cal easy magnetization axes are alternately arranged on the 
zigzag Mn-0 chain (upper figure), which can cause deviation 
of each Mn-spin direction from the b axis towards the local 
easy axis, resulting in the cycloidal deformation (lower fig- 
ure), (c) On the zigzag Mn-0 chain, the c-axis components 
of the Dzyaloshinskii-Moriya vectors are arranged in a stag- 
gered way, with which the Mn spins cant and rotate in the ab 
plane. Here (®) denotes the positive (negative) c-axis com- 
ponent of the Dzyaloshinskii-Moriya vector. Note that both 
the single-ion anisotropy with alternate easy magneti- 
zation axes and the Dzyaloshinskii-Moriya interaction Hdm 
give the equivalent spin canting or rotation as can been seen 
in lower figure of (b) and (c). 



mation within the ab plane. This is in contrast to what 
has been believed so far. Figure [7] depicts the real-space 
spin configuration of the -E-typc order calculated at T=0, 
which indeed shows an elliptically deformed afr-plane cy- 
cloid. 

There are two possible origins for this cycloidal de- 
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FIG. 9: (Color online) Temperature dependence of expecta- 
tion values for Hum an d for Jt—IA meV. At the tran- 
sition to the -E-type phase with lowering temperature, the 
value for "H^a abruptly decreases with a large jump, while 
that for Hdm slightly increases, indicating that the term H sia 
is a source of the cycloidal deformation. 



formation. One is the single- ion anisotropy Hg ia , and 
the other is the Dzyaloshinskii-Moriya interaction Hum ■ 
Let us consider a pure up- up- down- down spin structure 
shown in Fig. [SJa) in which the Mn spins are collincarly 
aligned parallel to the b axis, and examine what happens 
if we switch on H® m or Hum- The factor (— 1)*»>+ , » in 
Hsi& given in Eq. (6) implies that the local easy magne- 
tization axes are alternately arranged along the in-plane 
Mn-0 chain because of the staggered orbital ordering. 
Since the occupied e g orbital is d 3 ^2_ r 2 or d 3ri 2_ r 2 , the 
easy magnetization axis is along the or rji axis as shown 
in upper figure of Fig. [HJb) Therefore, in the presence of 
Tibial direction of each Mn spin deviates from the b axis 
toward £j or rji axis to form an clliptically deformed ab- 
plane cycloid as shown in lower figure of Fig. HJb) . On the 
other hand, in the presence of Hum, the c-axis compo- 
nents of the Dzyaloshinskii-Moriya vectors are arranged 
in a staggered way on the in-plane zigzag Mn-0 chains, 
which can also cause canting of the Mn spins to form a 
cycloidal rotation as shown in Fig. |8tc) . 

To identify the origin of the cycloidal deformation of 
the .E-type spin structure, we calculate temperature- 
dependence of the expectation values for Hum and H<f a 
at Jh=1.4 meV (sec Fig. [9]). Here H D b M represents the 
Dzyaloshinskii-Moriya interaction associated with the 
vectors d[j on the in-plane Mn-0 bonds, i.e., 



H 



ab 
DM 



^ ' ' {Si X Si+ X ) + ^ ' ' {Si X Si- 



(24) 

At the transition to the E-type phase with lowering tem- 
perature, the energy for H^ abruptly decreases with a 
large jump, while that for W^m increases slightly, indi- 
cating the single-ion anisotropy H^ ia or the alternation 
of the in-plane easy magnetization axes as its origin. 

The validity of the noncollincarity in the E-type phase 
predicted in the classical spin model is justified by the 
large Mn S—2 spins. The quantum fluctuation is almost 



suppressed in i?Mn03, and this is the reason why our 
model has successfully described a lot of experimental 
results for i?Mn03. More concretely, our model has been 
established by quantitative reproductions of the phase di- 



33.34.53 



and the electromagnon optical spectr; 



,16 



agrams 

These facts strongly support robustness of the predicted 
noncollinear E-type order. Furthermore, the predicted 
extent of the deformation sensitively depends on the 
strength of the single-ion anisotropy or the parameter 
value of E. Our choice of the parameter E=0.25 meV 
has quantitatively reproduced the area of the sinusoidal 
collinear phase in the phase diagram a 33 ' 34 ' 53 , the thresh- 
old magnetic field of the field-induced P reorientation^, 
the ellipticity of the cycloidal spin structures 3 ^, and the 
electromagnon spectral as experimentally observed, all 
of which are also sensitive to the value of E. These 
facts guarantee the validity of our parameter choice and, 
hence, the predicted degree of the noncollinear deforma- 
tion. 

There are already several neutron-scattering studies for 
i?Mn0 3 with R=Y. Ho, and Er^SI, but no apparent 
noncollinear deformation has been observed. One might 
think that the noncollinear deformation shown in Fig. 
is large enough to be detected experimentally. However 
it should be mentioned that most of the previous exper- 
iments failed to observe the real E-type phase. Indeed 
they reported an incommensurate wave numbers contra- 
dicting obviously to the E-type state. Thus far, only one 
experiment by Munoz et al. probably measured the real 
E-type phase with commensurate gh=0.57r in HoMnOy^, 
but it was performed for powder samples and without any 
electrical poling procedures. In fact, comparison between 
calculated Rietveld pattern for the predicted noncollinear 
E-type state and that for the pure collinear E-type state 
in the case of powder sample revealed that only slight dif- 
ferences appear in the intensity of peaks at (0 0.5 1) and 
(1 0.5 1) with Pbnm setting. The differences are ^10% 
for the former peak, while ~12% for the latter peak^. 
In turn, Munoz et al. compared their observed and cal- 
culated patterns, which also shows approximately 10% 
errors for both peaks (see Fig. 11(c) of Refi^). These er- 
rors suggest that accuracy of their measurement is not 
enough, or their analysis assuming the collinear E-type 
state is not appropriate. Now we would like to suggest 
that there is a possibility to achieve better agreement be- 
tween the observed and calculated patterns if they per- 
form the analysis assuming the predicted noncollinear 
E-type state. The expected signal of the noncollincar- 
ity is very small, and in order to detect such a small 
difference, a careful sample synthesis, measurement, and 
analysis are required. Quite recently, Ishiwata and his 
coworkers have succeeded in synthesizing the pcrovskite 
YM11O3 single crystals 3 ^, but their sizes (~0.5 mm) are 
not large enough for a neutron-scattering experiment. 

In the cycloidally deformed E-type spin order, both 
{S ■ S)-type and (S x S')-type mechanisms contribute 
to the ferroelectric P. Concerning the former mecha- 
nism, with alternate large and small spin rotation angles 
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FIG. 10: (Color online) In the spirally deformed up-up-down- 
down structure, both (S ■ S')-type and (S x S)-type mecha- 
nisms contribute to the ferroelectric P. Concerning the for- 
mer mechanism, alternate large and small spin turn angles 
cause a uniform shift of the O ions to strengthen and weaken 
the ferromagnetic exchanges through decreasing and increas- 
ing the Mn-O-Mn bond angle, respectively. Shifts of the O 
ions due to the (S ■ S)-type magnetostriction are shown by 
small (gray) arrows. Large (red and blue) arrows represent 
the dominant b-axis components of the Mn spins. 



or with dominant up- up- down- down spin 6-axis compo- 
nents, the O ions between nearly (anti)parallel Mn-spin 
pairs shift negatively (positively) to strengthen (weaken) 
the ferromagnetic exchanges through increasing (decreas- 
ing) the Mn-O-Mn bond angle, which results in the uni- 
form electric polarization. As shown in Fig. [TU1 the oxy- 
gen displacements on the zigzag Mn-0 chain along the 
x (y) axis contribute to the uniform electric polarization 
7r x (7r y ) pointing in the —y (x) direction, respectively. 
Consequently, the P-type spin order generates the fer- 
roelectric Pg parallel to the a axis as a sum of 7r x ||— y 
and n y \\+x. Wc also expect a small but finite (S x S) 
contribution Pas due to the cycloidal deformation. 

Now we discuss roles of the Peierls-type spin-phonon 
coupling and the weak ferromagnetic coupling J a along 
the a axis for realization of the P-typc spin order. The E- 
type order appears when the antiferromagnetic coupling 
Jb along the b axis is sufficiently strong. In order to 
specify the origin of the P-typc order, we consider the 
following two-dimensional classical Heisenberg model, 

i 

+ Jb Si ■ + J a Si ■ ffj+a, (25) 



with ferromagnetic Jy = J & + J' ab Si.j<0, antiferromag- 
netic Jb>0, and weakly ferromagnetic J o <0 [see also 
Fig. HIT a')]. This is a simplified model to understand 
the P-type order where the single-ion anisotropy and the 
Dzyaloshinskii-Moriya interaction are removed from the 
Hamiltonian When J' ab =0 and J a =0, three kinds 
of magnetic states are degenerate in the limit of strong 




/ / / / 
/ > / > 
> / > > 











t— % 












(b) (c) (d) 



JdJ_ 



Jab-0 
Ja=0 



(b) (C) 

Jab^O 
Ja=0 



(d) 



(c) 



(b) 



Jab'*0 
Ja*0 



FIG. 11: (Color online) (a) Spin exchanges of the classical 
Heisenberg model given by Eq. (|25[) . In the limit of large 
antiferromagnetic Ji,, three kinds of magnetic states are de- 
generate, i.e., (b) (pure) iJ-type, (c) stripe, and (d) 90° spiral 
states. The nearest neighbor bonds with energy gain (cost) of 
— \J a b\S 2 (+\J a b\S 2 ) are shown by thick (dotted) lines. The 
i5-type state becomes energetically stabilized by bond alter- 
nation shown in (b), while the stripe state by that shown in 
(c). (e) Energy diagram of these three states (see text). 



Jb, i.e., the (pure) P-type, stripe, and 90° spiral states, 
which are shown in Figs. [rTT bVfd). respectively. In all 
these three states, the spins align antiferromagnetically 
along the b axis. In the P-type and the stripe states, 
there are energy gains (costs) of —\J a b\S 2 {+\J a b\S 2 ) 
associated with the nearest neighbor ferromagnetic ex- 
changes J a b on the bonds connecting the parallel (an- 
tiparallel) spin pairs as indicated by thick (dotted) lines. 
These gains and costs are perfectly canceled out in total 
for both cases. On the other hand, there is neither gain 
nor cost in energy associated with J a b in the 90° spiral 
state because Si ■ Si is always zero in this state. Namely 
when only the spin-exchange interactions J a b and Jb are 
considered, the energies of these three states are all iden- 
tical resulting in their degeneracy. Then if we incorporate 
the spin-lattice coupling by taking finite J' ab , a simul- 
taneous bond alternation sets in to lift the degeneracy 
by modulating the nearest neighbor ferromagnetic ex- 
changes. The E-type and the stripe states become lower 
in energy by the bond alternations shown in Figs. Illf b) 
and (c), respectively. Here the ferromagnetic exchanges 
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FIG. 12: (Color online) (a) Calculated correlation function 
J 7 y (fc, T) given by Eq. (f26|) for the -E-type state at J&=2.4 
meV and fcBT=0.5 meV, which shows sharp peaks at fe=(±7T, 
±7r, 0). Magnitudes of these peaks are all l.OxlO -4 meV 2 . 
(b) Modulations of the nearest-neighbor ferromagnetic ex- 
changes. The ferromagnetic coupling is strengthened (weak- 
ened) by AJ a (,~0.01 meV on the thick (thin) bonds. 

are strengthened (weakened) on the thick (dotted) bonds. 
In contrast, the energy of the 90° spiral state does not 
change by bond alternations. 

The remaining degeneracy of the £7- type and the stripe 
states is eventually lifted by the weak ferromagnetic in- 
teraction J a . Namely the E-type state with ferromagnct- 
ically aligned spins along the a axis becomes stabilized by 
J a , while the stripe state with the staggered spin align- 
ment along the a axis does not. As we have discussed, all 
of the three terms in the Hamiltonian (|2"5|) arc indispens- 
able for the -E-type order, and other terms in the original 
Hamiltonian i.e., the single-ion anisotropy and the 
Dzyaloshinskii-Moriya interaction, rather favor the non- 
collinear spin alignment. Therefore we conclude that the 
£?-type order in i?Mn03 is a consequence of the three in- 
gredients, i.e., the large next neighbor antiferromagnetic 
coupling Jh, the Peierls-type spin-lattice coupling J' a (,Sij, 
and the weak ferromagnetic coupling J a . 

To confirm that the bond alternation shown in 
Fig. UTT b) is actually realized in the E-type phase, we 
calculate the following correlation functions at Jb=2A 
meV and /cbT'=0.5 meV, 

J 77 ,(fc,T) = ^J2(AJ iii+1 AJ jtj+Y )e ik < r ^\ (26) 

i,i 

for (7, j')=(x, x), (y, y) and (x, y). Here AJ t j=J' ab 5ij. 
We find that all these correlation functions are identi- 
cal. They have sharp peaks at fe=(±7r, ±7r, 0), and 
their magnitudes are all equal to l.OxlO -4 meV 2 as 
shown in Fig. [T2Ta). This means that the bond alter- 
nation with J a f, ± A J a b is indeed realized where A J ao ~ 

J J 77 '(7r, 7r, 0) is approximately 0.01 meV, i.e., 1.25 % of 
the original value |J a (,|=0.8 meV [see Fig. IT2Tb)]. 

C. Spiral Spin Phases 

In Figs. rHfl' a) and (b), we show calculated temperature 
profiles of specific heat C S (T), spin helicity /i 7 (T), and 



polarization P 7 (T) due to the (S ■ S')-type magnetostric- 
tion for (a) J b =0.7 meV and (b) J 6 =1.2 meV. For J b =0.7 
meV, three successive phase transitions take place with 
lowering temperature. From the paramagnetic phase, 
the system first enters into the sinusoidal collincar phase. 
Subsequently the system enters into the fee-plane spiral 
phase, and finally into the afe-plane spiral phase. On 
the other hand, when J =l.2 meV, the system exhibits 
only two phase transitions among three phases, i.e., the 
paramagnetic, the sinusoidal collinear, and the fee-plane 
spiral phases, whereas the afe-plane spiral phase does not 
appear. 

Interestingly, we see that P a (T) in the afe-plane spiral 
phase for J b =0.7 meV is extrapolated to ^500 fiC/m 2 
at T->0 [lower panel of Fig. [THTa)]. while to zero in 
the fee-plane spiral phase for J;, =1.2 meV [lower panel 
of Fig. [I3](b)]. This indicates a finite contribution to 
the ferroelectric polarization from the (S ■ S)-type mag- 
netostriction in the afe-plane spiral phase. Moreover we 
find that this (S ■ S) contribution can be comparable to 
or even larger than the (S x S) contribution as discussed 
later. This is surprising because only the (S x S^-type 
magnetostriction has been considered as an origin of the 
ferroelectric order in the spiral spin phase thus far. Con- 
trastingly the ferroelectric order in the fee-plane spiral 
phase is purely of (S x S) origin with no (S ■ S) contri- 
bution. 

This (S ■ S) contribution can explain puzzling experi- 
mental results. It has been known that the amplitude of 
the ferroelectric P in the a&-plane spiral phase is much 
larger than that in the frc-plane spiral phase. For in- 
stance, P in the a6-plane spiral phase of DyMnC>3 under 
H\\b is 2.5 times larger than P in the &c-plane spiral 
phase at £f=0l&2i. Moreover, in Euo. 6 Yo.4Mn0 3 , the 
P || a at H=0 is approximately ten times larger than P\\c 
under H\\a^. Importantly the latter example excludes 
the possible influence of /-electron moments as its origin 
because of their absence in Eu 3+ and Y 3+ ions, which 
enables us to highlight the roles of spin-phonon coupling 
on the polarization behavior. These observations have 
been a puzzle since we expect nearly identical strength 
of the (S x S)-type magnetostriction in the a6-plane and 
fee-plane spiral phases. 

The (S ■ S) contribution in the afe-plane spiral phase 
can be understood by a combined function of the 
Dzyaloshinskii-Moriya interaction and the symmetric 
magnetostriction (see also Fig. IT4|) . On the in-plane 
zigzag Mn-0 chains, the c-axis components of the 
Dzyaloshinskii-Moriya vectors are arranged in the stag- 
gered way. Under this circumstance, the spin rotation 
angles in the a&-plane spiral phase become subject to an 
alternate modulation. Then the O ions between two spins 
with a smaller angle of (f)—A<p (a larger angle of 4>+A<j>) 
shift negatively (positively) to strengthen (weaken) the 
ferromagnetic exchange through increasing (decreasing) 
the Mn-O-Mn bond angle. These shifts generate a uni- 
form component resulting in the ferroelectric polariza- 
tion. In fact, the spin rotation angles in the frc-planc 




FIG. 13: (Color online) (a) Temperature profiles of specific heat C S (T), spin helicity /i 7 (T), and polarizations due to the 
(S ■ S)-type magnetostriction P-,(T) for Jb— 0.7 meV. (b) Those for Jt=1.2 meV. 




FIG. 14: Alternation of the spin turn angles in the a&-plane 
spiral state due to the staggered Dzyaloshinskii-Moriya vec- 
tors is depicted in an exaggerated manner where (®) de- 
notes the positive (negative) c-axis component of the vector. 
Induced shifts of the O ions due to the (S ■ S)-type magne- 
tostriction are shown by gray arrows. 



spiral phase are also subject to the alternate modula- 
tion because of the staggered a-axis components of the 
Dzyaloshinskii-Moriya vectors. However the induced O 
shifts are opposite between neighboring ab planes, which 
results in their perfect cancellation. 

The (S ■ S) mechanism in the afr-plane spiral phase 
is triggered by the alternate spin-angle modulation due 
to the staggered Dzyaloshinskii-Moriya vectors, and thus 
is a higher-order effect of the Dzyaloshinskii-Moriya in- 
teraction. One might think that this (S ■ S) contribu- 
tion Pg cannot be larger than the (S x S) contribu- 
tion Pas because the latter is a direct consequence of 
the Dzyaloshinskii-Moriya coupling. But this is not cor- 
rect, and the Ps can be much larger than the Pas m 
reality. This is because the inherent "large" staggered 
Dzyaloshinskii-Moriya vectors cause the "large" alternate 
spin-angle modulation, which results in the generation of 
large Pas through the symmetric (S • S')-type magne- 
tostriction. On the other hand, the spiral spin order 
causes weak "fcrri-components" of the Dzyaloshinskii- 
Moriya vectors via the antisymmetric (S x S^-type mag- 
netostriction through inducing uniform oxygen shifts or 



uniform deformations of the electron distribution. How- 
ever, these (S x S)-induced ferri-components are so small 
that Pas tends to be small. 

Recently Malashevich and Vanderbilt found in their 
full first-principles calculation that the ferroclcctricity in 
TbMn03 (&c-plane spiral) cannot be explained by the 
simple (S x S) mechanism onl y 69 ' 70 , and they suggested 
the presence of other coupling. However we should note 
that the (S ■ S) magnetostriction proposed in the present 
paper cannot be "other coupling" suggested by them be- 
cause the (S ■ S) mechanism works only in the a6-planc 
spiral phase but not in the 6c-plane spiral phase. The 
ferroelectricity in TbMn03 still contains some puzzles in 
its mechanism, which should be uncovered in the future 
study. 

We also calculate the J^-dependence of the (S x S) 
contribution Pas at T— >0 within the spin-current model 
where Pas is proportional to the spin helicity h = 
J2<ij> Si x Sj. Here we use the temperature profile 
of the spin helicity h 7 (T) in Eq. (17), which is obtained 
by the Monte-Carlo simulation. Note that the Pas is 
not calculated from the oxygen positions Si j in contrast 
to the (S ■ S) contribution Ps. This is because there 
are two contributions to PaSj i.e, the electronic and the 
lattice-mediated contributions as pointed out by Mala- 
shevich and Vanderbil t 69 ' 70 . Since the electronic contri- 
bution (contribution from the deformed electron-clouds 
around atoms) cannot be evaluated in our spin model, we 
calculate the value of Pas at T— >0 from h 7 (T— >0) shown 
in Fig. [15] using the proportional relation between them. 
The proportionality factor can be evaluated by compar- 
ing the calculated value of h a (T— >0) and the experimen- 
tal value of P in the 6c-plane spiral phase because the ob- 
served P in the 6c-plane spiral phase is purely of (S x S) 
origin. The calculated value of h a (T— H)) is nearly con- 
stant in the &c-plane spiral phase and is 0.875 at J b =1.2 
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FIG. 15: (Color online) Calculated Jb dependence of the mag- 
netic wave number qi, and the 7-axis components of the spin 
helicity ft 7 (l~ a , b, c) at T— 5>0. 
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FIG. 16: (Color online) Jt, dependence of polarizations at 
T— K), i.e., (S ■ S) contribution Ps, (S x 5) contribution 
Pas, and experimentally measured P in Eui-a; Y^MnOa and 
Yi-yLuyMnOg^. The summation Ps+Pas reproduces the 
experimental P well. 



meV, while the observed Pas(P— >0) is ~455 /iC/m 2 for 
Eui^Y^MnOs with x=0.75^. These gives the propor- 
tionality factor of 520 /iC/m 2 . 

In Fig. [16l we plot the calculated J^-dcpcndence of 
-Ps, -Pas, an d Ps + -Pas- We also plot experimentally 
measured P for the solid solutions Eui-^Y^MnOa and 
Yi_ a Lu y Mn03 for comparison^, whose P originates 
purely from the Mn-spin order because of the absence 
of / moments. Effective r# and J b of these solid solu- 
tions are evaluated by interpolations. We find that the 
sum Ps+Pas reproduces well the experimental P. In 
particular, our calculation gives constant P in the P- 
type phase in agreement with the experiment. It should 
be emphasized that only the elastic constant K is an un- 
known parameter in our model, and once we determine 
its value so as to reproduce the experimental P in the 
P-type phase, the behaviors of P in the spiral phases 
are reproduced almost perfectly. Moreover it turns out 
that the (S ■ S) contribution Ps can be comparable to 



or even larger than the (S x S) contribution Pas in the 
a6-plane spiral phase. This explains why the P in the 
afe-plane spiral phase is much larger than that in the fee- 
plane spiral phase. The 2.5 times larger P\\a under H\\b 
than P\\c at H=0 in DyMn0 3 is ascribed to this (S ■ S) 
contribution. We expect that the (S ■ S) contribution 
is 1.5 times larger than the (S x S) contribution in the 
afe-planc spiral phase of DyMnC>3. 



D. Coexisting States 

Next we discuss a certain kind of metastable incom- 
mensurate spin state and its possible coexistence with the 
commensurate P-type state in the large Jb region. The 
spin-phonon coupling or the (S ■ S^-type magnetostric- 
tion make the transition between the incommensurate 
spiral and the P-type phases of strong first order. Con- 
sequently, some incommensurate states have deep local 
energy minima even in the P-type phase although the 
energy comparison gives the transition line as indicated 
by the solid line in Fig. [5] This can result in the realiza- 
tion of a metastable incommensurate spin state trapped 
in a local energy minimum or its coexistence with the 
P-type state. They can easily occur in reality since the 
system enters into the P-type phase necessarily via the 
incommensurate sinusoidal collinear phase with lowering 
temperature. 

For 1.4<J& (meV)<2.5, since both the P-type and 
the incommensurate states have energy minima, even 
the replica-exchange Monte-Carlo calculation sometimes 
fails to reach the real lowest-energy state. Thus we per- 
form calculations starting with certain initial spin con- 
figurations. Note that all other calculations are started 
with random spin configurations. We chose (A) P-type 
spin configuration obtained for p,=2.4 meV and (B) in- 
commensurate spiral states obtained by switching off the 
spin-lattice coupling or by setting J' ab =0 as initial con- 
figurations. In both cases, we perform the Monte-Carlo 
sampling after sufficient thermalization. 

In Fig. 1171 we show the real-space spin configuration 
obtained in the calculation starting with (B). We can see 
that small incommensurate (<#,=(). 4587r) spiral regimes 
exist in the background commensurate (<#,=(). 57r) phase. 
Interestingly we find that these incommensurate regimes 
emerge periodically to form a stripe structure, indicating 
possible realization of magnetic discommensulation. 

In fact, there exists an apparent contradiction in the 
neutron-diffraction results for PMnC>3 with small R 
ions^— . Most of the previous experiments reported 
incommensurate wave numbers g&^0.437r in HoMnOg^, 
YMnOg^, and ErMnO^. Moreover one of the reports 
claimed that the magnetic structure in YMnC>3 is simul- 
taneously incommensurate and collinear even down to 
the lowest temperature of 1.7 K— . On the other hand, 
one experiment reported a commensurate wave number 
of ^b=0.57r in HoMnOg^. This puzzle can be solved by 
considering the presence of above metastable incommen- 
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suratc spin state. We calculate the spin-correlation func- 
tions in the momentum space for the above metastablc 
solution. We find that only the spin &-axis component has 
sharp peaks at <7t,=±0.4587r, while the other components 
have very small peaks as shown in Fig. [TST a) . This seems 
as if the spin structure were incommensurate collinear. 
For comparison we also display the calculated spin cor- 
relation function for the commensurate S-type state in 
Fig. ll8f b) where the spin fe-axis (a-axis) components have 
large (small) peaks at qb=±0.5ir. Accordingly the ob- 
served incommensurate wave numbers and the claimed 
incommensurate collinear state in YM11O3 can be at- 
tributed to the metastable incommensurate state, while 
a report of the commensurate qb=0.57r in HoMnO^ can 
be ascribed to the pure E-type state. The metastable 
incommensurate spin state and its coexistence with the 
-E-type state should be seriously considered also when we 
interpret the experimental results for .RMJ1O3 with R=Y. 
Ho, ...,Lu, such as strange electromagnon spectra in the 
THz optical spectroscopy-^. 



10 20 30 40 50 60 70 80 90 
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FIG. 17: (Color online) Real-space spin configuration ob- 
tained in the Monte- Carlo simulation starting with an in- 
commensurate spiral spin configuration as the initial state at 
J b =2A meV (see text). The spin 7-axis components (7=0, b, 
c) are plotted for two kinds of ab plane, i.e., Zi~0 and Zi—0.5. 




FIG. 18: (Color online) (a) Calculated spin-correlation func- 
tions in the momentum space at Jt=2.4 meV for a metastable 
incommensurate state obtained in the Monte-Carlo simu- 
lation starting with the incommensurate spiral state (see 
text). Their peaks are located at gi,=±0.4587r. (b) Those 
for the commensurate iJ-type state whose peaks are located 
at qb—±0.5n. Here 5* 7 (fc) denotes the correlation function for 
the spin 7-axis components given by Eq. (|19|l . 



V. SUMMARY 



In summary, we have theoretically studied the origins 
and nature of the multiferroic phases and the magne- 
toelectric coupling in i?Mn03 by using a realistic spin 
model including the Peierls-type spin-phonon coupling, 
which can successfully reproduces the entire phase dia- 
gram of i?Mn03. We have revealed the cooperative con- 
tributions of symmetric (S ■ S^-type and antisymmetric 
(S x 5)-typc magnetostrictions to the ferroelectricity in 
the afr-plane spiral phase. This large (S ■ S) contribu- 
tion to ferroelectricity is expected and should be seri- 
ously considered also in other spin-spiral-based multifer- 
roics, for which the (S x S^-type magnetostriction has 
been believed to be a unique origin of the ferroelectric 
polarization. We have also uncovered the cycloidal spin 
deformation in the -E-type phase due to the alternate ar- 
rangement of the easy-magnetization axes on the in-plane 
zigzag MnO chain. We expect that this cycloidal defor- 
mation would be detected in a future neutron-scattering 
experiment. A metastable incommensurate spin state 
and its coexistence with the E-type state have been also 
found on the verge of the phase boundary between the 
E-type and spiral states. On these basis, a puzzle in 
the neutron-scattering experiments for i?=Y, Ho, and Er 
have been resolved. Our model gives a firm basis for 
studying and controlling the intriguing cross-correlation 
phenomena in /2M11O3. Moreover the crucial roles of the 
spin-phonon coupling we have demonstrated by taking 
.RM11O3 are not specific to this manganite system, but 
are relevant to all of the multiferroic materials. 



14 



Acknowledgment 

The authors are grateful to Y. Tokura, S. Ishiwata, 
F. Kagawa, D. Okuyama, and T. Arima for discussions. 
This work was supported by Grant-in-Aid for Scien- 
tific Research (Grants No. 22740214, No. 21244053, 
No. 17105002, No. 19048015, and No. 19048008), 



Global-COE Program ("Physical Sciences Frontier") and 
NAREGI Project from the Ministry of Education, Cul- 
ture, Sports, Science and Technology of Japan, and by 
Funding Program for World-Leading Innovative R&D on 
Science and Technology (FIRST Program) on "Quantum 
Science on Strong Correlation" . 



1 I. E. Dzyaloshinskii, Zh. Eksp. Teor. Fiz. 33, 881 (1959) 
[Sov. Phys. JETP 10, 628 (I960)]. 

2 D. N. Astrov, Zh. Eksp. Teor. Fiz. 38, 984 (1960) [Sov. 
Phys. JETP 11, 708 (I960)]. 

3 Magnetoelectric Interaction Phenomena in Crystals, edited 
by A. J. Freeman and H. Schmid (Gordon and Breach, Lon- 
don, 1975). 

4 T. Kimura, T. Goto, H. Shintani, K. Ishizaka, T. Arima, 
and Y. Tokura, Nature (London) 426, 55 (2003). 

5 H. Schmid, Ferroelectrics 162, 317 (1994). 

6 N. A. Hill, J. Phys. Chem. B 104, 6694 (2000). 

7 M. Fiebig, J. Phys. D: Appl. Phys. 38, R123 (2005). 

8 D. I. Khomskii, J. Magn. Magn. Mater. 306, 1 (2006). 

9 Y. Tokura, Science 312, 1481 (2006). 

10 W. Eerenstein, N. D. Mathur, and J. F. Scott, Nature 
(London) 442, 759 (2006). 

11 S.-W. Cheong and M. Mostovoy, Nat. Mater. 6, 13 (2007). 

12 Y. Tokura, J. Magn. Magn. Mater. 310, 1145 (2007). 

13 A. Pimenov, A. A. Mukhin, V. Yu. Ivanov, V. D. Travkin, 
A. M. Balbashov, and A. Loidl, Nat. Phys. 2, 97 (2006). 

14 N. Kida, Y. Takahashi, J. S. Lee, R. Shimano, Y. Ya- 
masaki, Y. Kaneko, S. Miyahara, N. Furukawa, T. Arima, 
and Y. Tokura, J. Opt. Soc. Am. B 26, A35 (2009). 

15 R. ValdesAguilar, M. Mostovoy, A. B. Sushkov, C. L. 
Zhang, Y. J. Choi, S-W. Cheong, and H. D. Drew, Phys. 
Rev. Lett. 102, 047203 (2009). 

16 M. Mochizuki, N. Furukawa, and N. Nagaosa, Phys. Rev. 
Lett. 104, 177206 (2010). 

17 M. Mochizuki and N. Nagaosa, Phys. Rev. Lett. 105, 
147202 (2010). 

18 T. Kimura, G. Lawes, T. Goto, Y. Tokura, and A. P. 
Ramirez, Phys. Rev. B 71, 224425 (2005). 

19 M. Tokunaga, Y. Yamasaki, Y. Onose, M. Mochizuki, N. 
Furukawa, and Y. Tokura, Phys. Rev. Lett. 103, 187202 
(2009). 

20 N. Abe, K. Taniguchi, S. Ohtani, T. Takenobu, Y. Iwasa, 
and T. Arima, Phys. Rev. Lett. 99, 227206 (2007). 

21 H. Murakawa, Y. Onose, F. Kagawa, S. Ishiwata, Y. 
Kaneko, and Y. Tokura, Phys. Rev. Lett. 101, 197207 
(2008). 

22 M. Mochizuki and N. Furukawa, Phys. Rev. Lett. 105, 
187601 (2010). 

23 T. Goto, T. Kimura, G. Lawes, A. P. Ramirez, and Y. 
Tokura, Phys. Rev. Lett. 92, 257201 (2004). 

24 F. Kagawa, M. Mochizuki, Y. Onose, H. Murakawa, 
Y. Kaneko, N. Furukawa, and Y. Tokura, Phys. Rev. 
Lett. 102, 057604 (2009). 

25 F. Schrettle, P. Lunkenheimer, J. Hemberger, V. Yu. 
Ivanov, A. A. Mukhin, A. M. Balbashov, and A. Loidl, 
Phys. Rev. Lett. 102, 207208 (2009). 

26 H. Katsura, N. Nagaosa, and A. V. Balatsky, Phys. Rev. 
Lett. 95, 057205 (2005). 



27 I. A. Sergienko and E. Dagotto, Phys. Rev. B 73, 094434 
(2006). 

28 M. Mostovoy, Phys. Rev. Lett. 96, 067601 (2006). 

29 M. Kenzelmann, A. B. Harris, S. Jonas, C. Broholm, J. 
Schefer, S. B. Kim, C. L. Zhang, S.-W. Cheong, O. P. Vajk, 
and J. W. Lynn, Phys. Rev. Lett. 95, 087206 (2005). 

30 Y. Yamasaki, H. Sagayama, T. Goto, M. Matsuura, K. 
Hirota, T. Arima, and Y. Tokura, Phys. Rev. Lett. 98, 
147204 (2007). 

31 T. Arima, A. Tokunaga, T. Goto, H. Kimura, Y. Noda, 
and Y. Tokura, Phys. Rev. Lett. 96, 097202 (2006). 

32 Y. Yamasaki, H. Sagayama, N. Abe, T. Arima, K. Sasai, 
M. Matsuura, K. Hirota, D. Okuyama, Y. Noda, and Y. 
Tokura, Phys. Rev. Lett. 101, 097204 (2008). 

33 M. Mochizuki, and N. Furukawa, J. Phys. Soc. Jpn. 78, 
053704 (2009). 

34 M. Mochizuki, and N. Furukawa, Phys. Rev. B 80, 134416 
(2009). 

35 I. A. Sergienko, C. Sen, and E. Dagotto, Phys. Rev. Lett. 
97, 227204 (2006). 

36 S. Picozzi, K. Yamauchi, B. Sanyal, I. A. Sergienko, and 
E. Dagotto, Phys. Rev. Lett. 99, 227201 (2007). 

37 K. Yamauchi, F. Freimuth, S. Blugel, and S. Picozzi, Phys. 
Rev. B 78, 014403 (2008). 

38 S. Ishiwata, Y. Kaneko, Y. Tokunaga, Y. Taguchi, T. 
Arima, and Y. Tokura, Phys. Rev. B 81, 100411(R) (2010). 

39 S. Ishiwata, Y. Tokunaga, Y. Taguchi, and Y. Tokura, J. 
Am. Chem. Soc. 133, 13818 (2011). 

40 B. Lorenz, Y. Q Wang, and C. W. Chu, Phys. Rev. B 76, 
104405 (2007). 

41 V. Yu. Pomjakushin, M. Kenzelmann, A. Donni, A. B. 
Harris, T. Nakajima, S. Mitsuda, M. Tachibana, L. Keller, 
J. Mesot, H. Kitazawa, and E. Takayama-Muromachi, New 
J. of Phys. 11, 043019 (2009). 

42 T. A. Kaplan and S. D. Mahanti, Phys. Rev. B 83, 174432 
(2011). 

43 T. Goto, Y. Yamasaki, H. Watanabe, T. Kimura, and Y. 
Tokura, Phys. Rev. B 72, 220403(R) (2005). 

44 J. Hemberger, F. Schrettle, A. Pimenov, P. Lunkenheimer, 
V. Y. Ivanov, A. A. Mukhin, A. M. Balbashov, and A. 
Loidl, Phys. Rev. B 75, 035118 (2007). 

45 Y. Yamasaki, S. Miyasaka, T. Goto, H. Sagayama, T. 
Arima, and Y. Tokura, Phys. Rev. B 76, 184418 (2007). 

46 V. Yu. Ivanov, A. A. Mukhin, V. D. Travkin, A. S. 
Prokhorov, A. M. Kadomtsev, Yu. F. Popov, G. P. 
Vorobev, K. I. Kamilov, and A. M. Balbashov, J. Magn. 
Magn. Mater. 300, el30 (2006). 

47 V. Yu. Ivanov, A. A. Mukhin, V. D. Travkin, A. S. 
Prokhorov, Yu. F. Popov, A. M. Kadomtseva, G. P. 
Vorobev, K. I. Kamilov, A. M. Balbashov, Phys. Status 
Solidi B 243, 107 (2006). 

48 T. Kimura, S. Ishihara, H. Shintani, T. Arima, K. T. 



15 



Takahashi, K. Ishizaka, and Y. Tokura, Phys. Rev. B 68, 
060403(R) (2003). 

49 T. Hotta, M. Moraghebi, A. Feiguin, A. Moreo, S. Yunoki, 
and E. Dagotto, Phys. Rev. Lett. 90, 247203 (2003); S. 
Dong, R. Yu, S. Yunoki, J.-M. Liu, and E. Dagotto, Phys. 
Rev. B 78, 155121 (2008). 

50 T. A. Kaplan, Phys. Rev. B 80, 012407 (2009). 

51 L. X. Hayden, T. A. Kaplan, and S. D. Mahanti, Phys. 
Rev. Lett. 105, 047203 (2010). 

52 N. Furukawa, and M. Mochizuki, J. Phys. Soc. Jpn. 79, 
033708 (2010). 

53 M. Mochizuki, N. Furukawa, and N. Nagaosa, Phys. Rev. 
Lett. 105, 037205 (2010). 

54 J. A. Alonso, M. J. Martmez-Lope, M. T. Casais, and M. 
T. Fernandez-Diaz, Inorg. Chem. 39, 917 (2000). 

55 I. Dzyaloshinsky, J. Phys. Chem. Solids 4, 241 (1958). 

56 T. Moriya, Phys. Rev. Lett. 4, 228 (1960). 

57 T. Moriya, Phys. Rev. 120, 91 (1960). 

58 I. Solovyev, N. Hamada, and K. Terakura, Phys. Rev. Lett. 
76, 4825 (1996). 

59 S. Picozzi, K. Yamauchi, G. Bihlmayer, and S. Blugcl, 
Phys. Rev. B 74, 094402 (2006). 

60 K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. 65, 1604 
(1996). 

61 Y. Tokura, Reports on Progress in Physics 69, 797 (2006). 

62 Y. Tomioka and Y. Tokura, Phys. Rev. B 70, 014432 



(2004). 

Successive emergence of these three magnetic phases was 
also observed in other magnetoelectric system, CuFeC>2; T. 
Kimura, J. C. Lashley, and A. P. Ramirez, Phys. Rev. B 
73, 220401 (2006). 

A. Munoz, M. T. Casais, J. A. Alonso, M. J. Martmez- 
Lope, J. L. Martinez, and M. T. Fernandez-Diaz, Inorg. 
Chem. 40, 1020 (2001). 

H. W. Brinks, J. Rodriguez-Carvajal, H. Fjellvag, A. Kjek- 
shus, and B. C. Hauback, Phys. Rev. B 63, 094411 (2001). 
A. Munoz, J. A. Alonso, M. T. Casais, M. J. Martmez- 
Lope, J. L. Martinez, and M. T. Fernandez-Diaz, J. Phys.: 
Condens. Matter 14, 3285 (2002). 

F. Ye, B. Lorenz, Q. Huang, Y. Q. Wang, Y. Y. Sun, C. 

W. Chu, J. A. Fernandez-Baca, Pengcheng Dai, and H. A. 

Mook, Phys. Rev. B 76, 060402(R) (2007). 

T. Arima and D. Okuyama, private communication. 

A. Malashevich and D. Vanderbilt, Phys. Rev. Lett. 101, 

037210 (2008). 

A. Malashevich and D. Vanderbilt, Phys. Rev. B 80, 
224407 (2009). 

Y. Takahashi, S. Ishiwata, S. Miyahara, Y. Kaneko, N. 
Furukawa, Y. Taguchi, R. Shimano, and Y. Tokura, Phys. 
Rev. B 81, 100413 (2010). 



